Alpha diversity
19 August 2024
1 Setup
1.1 Libs
## Loading required package: carData## Loading required package: lattice## Loading required package: plyr## Loading required package: ggplot2##
## Attaching package: 'ggpubr'## The following object is masked from 'package:plyr':
##
## mutate##
## Attaching package: 'cowplot'## The following object is masked from 'package:ggpubr':
##
## get_legend## Warning in checkDepPackageVersion(dep_pkg = "TMB"): Package version inconsistency detected.
## glmmTMB was built with TMB version 1.9.11
## Current TMB version is 1.9.14
## Please re-install glmmTMB from source or restore original 'TMB' package (see '?reinstalling' for more information)## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'## Loading required package: stats4## This is DHARMa 0.4.6. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')1.2 Data
setwd("~/nicolagk@hawaii.edu - Google Drive/My Drive/Mosquito_business/mosquito_microbes/04.microbiome_analysis/04a.diversity")
ps.clean <- readRDS("../../02.process_asvs/ps.clean.trim.less.rds")
ps.clean ## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 161 taxa and 458 samples ]
## sample_data() Sample Data: [ 458 samples by 24 sample variables ]
## tax_table() Taxonomy Table: [ 161 taxa by 11 taxonomic ranks ]2 Alpha diversity calculations
df.all <- data.frame(estimate_richness(ps.clean, split=TRUE, measures=c("Shannon","InvSimpson","Observed")))
df.all$sample_name <- rownames(df.all)
df.all.div <- merge(df.all,samdf,by="sample_name") #add sample data
#shannon diversity divided by species richness
df.all.div$even <- df.all.div$Shannon/(log(df.all.div$Observed))
df.div <- df.all.div
#adding most recent sample counts
sums <- data.frame(sample_sums(ps.clean))
colnames(sums) <- c("lib_size_trim")
sums$sample_name <- row.names(sums)
df.div <- merge(df.div,sums,by="sample_name")
#write.csv(df.div,"div.trim.csv")3 Analyzing alpha metrics
3.1 All data
Very strange/interesting pattern where the water types that shouldn’t have a lot of diversity have a lot of OTUs, but not a lot of reads…
ggplot(df.div,aes(x=type,y=lib_size_trim,color=infusion))+
geom_boxplot()+
geom_jitter(position=position_jitterdodge())
3.2 Plotz
#just experimental types
df.div.exp <- subset(df.div,type=="Microbial Water"|type=="A.albopictus")
df.div.exp$inf.temp <- paste0(df.div.exp$infusion,df.div.exp$temperature)
df.div.exp$type <- sub("A.albopictus","Mosquitoes",df.div.exp$type)
df.div.exp$type <- sub("Microbial Water","Meso. water",df.div.exp$type)
df.div.exp$type <- factor(df.div.exp$type,levels=c("Mosquitoes","Meso. water"))3.2.1 Richness
df.div.exp.se <- summarySE(df.div.exp,measurevar="Observed",groupvars=c("inf.temp","infusion","temperature","type"))
##dots & error bars
gg.rich <- ggplot(df.div.exp,aes(x=infusion,y=Observed,fill=inf.temp,shape=inf.temp))+
geom_jitter(alpha=0.5,position=position_jitterdodge(jitter.width=0.2),color="gray")+
geom_errorbar(data=df.div.exp.se,aes(ymax=Observed+se,ymin=Observed-se),width=0.2,color="black",position=position_dodge(width=0.6))+
geom_point(data=df.div.exp.se,size=2.5,position=position_dodge(width=0.6))+
#geom_boxplot(outlier.shape=NA,alpha=0.5)+
ylim(4,34)+
facet_wrap(~type)+
theme_cowplot()+
ylab("ASV richness")+
xlab("Infusion")+
scale_x_discrete(labels=c("OL","SG","PW"))+
scale_fill_manual(values=c("#B5D69F","#7DBA54","#C488A2","#992559","#B6F6EE","#60eede"),name="Infusion, temp.",labels=c("OL, cool","OL, warm","SG, cool","SG, warm","PW, cool","PW, warm"))+
scale_shape_manual(values=c(22,24,22,24,22,24),name="Infusion, temp.",labels=c("OL, cool","OL, warm","SG, cool","SG, warm","PW, cool","PW, warm"))
gg.rich
##water by time
df.div.water <- subset(df.div.exp,type=="Meso. water")
gg.mw.rich.time <- ggplot(df.div.water,aes(x=infusion,y=Observed,fill=inf.temp))+
geom_boxplot()+
#geom_jitter(alpha=0.5,position=position_jitterdodge(jitter.width=0.2),color="gray")+
#geom_errorbar(data=df.div.exp.se,aes(ymax=Observed+se,ymin=Observed-se),width=0.2,color="black",position=position_dodge(width=0.6))+
#geom_point(data=df.div.exp.se,size=2.5,position=position_dodge(width=0.6))+
#geom_boxplot(outlier.shape=NA,alpha=0.5)+
facet_wrap(~day)+
theme_cowplot()+
ylab("ASV richness")+
xlab("Infusion")+
scale_x_discrete(labels=c("OL","SG","PW"))+
scale_fill_manual(values=c("#B5D69F","#7DBA54","#C488A2","#992559","#B6F6EE","#60eede"),name="Infusion, temp.",labels=c("OL_cool","OL_warm","SG_cool","SG_warm","PW_cool","PW_warm"))+
ggtitle("Mesocosm water")
gg.mw.rich.time
3.2.2 Simpson’s
df.div.exp.se <- summarySE(df.div.exp,measurevar="InvSimpson",groupvars=c("inf.temp","infusion","temperature","type"))
##dots & error bars
gg.simp <- ggplot(df.div.exp,aes(x=infusion,y=InvSimpson,fill=inf.temp,shape=inf.temp))+
geom_jitter(alpha=0.5,position=position_jitterdodge(jitter.width=0.2),color="gray")+
geom_errorbar(data=df.div.exp.se,aes(ymax=InvSimpson+se,ymin=InvSimpson-se),width=0.2,color="black",position=position_dodge(width=0.6))+
geom_point(data=df.div.exp.se,size=2.5,position=position_dodge(width=0.6))+
#geom_boxplot(outlier.shape=NA,alpha=0.5)+
facet_wrap(~type)+
theme_cowplot()+
ylab("Simpson's (inv.)")+
xlab("Infusion")+
scale_x_discrete(labels=c("OL","SG","PW"))+
scale_fill_manual(values=c("#B5D69F","#7DBA54","#C488A2","#992559","#B6F6EE","#60eede"),name="Infusion, temp.",labels=c("OL, cool","OL, warm","SG, cool","SG, warm","PW, cool","PW, warm"))+
scale_shape_manual(values=c(22,24,22,24,22,24),name="Infusion, temp.",labels=c("OL, cool","OL, warm","SG, cool","SG, warm","PW, cool","PW, warm"))
gg.simp
gg.mw.simp.time <- ggplot(df.div.water,aes(x=infusion,y=InvSimpson,fill=inf.temp))+
geom_boxplot()+
#geom_jitter(alpha=0.5,position=position_jitterdodge(jitter.width=0.2),color="gray")+
#geom_errorbar(data=df.div.exp.se,aes(ymax=Observed+se,ymin=Observed-se),width=0.2,color="black",position=position_dodge(width=0.6))+
#geom_point(data=df.div.exp.se,size=2.5,position=position_dodge(width=0.6))+
#geom_boxplot(outlier.shape=NA,alpha=0.5)+
facet_wrap(~day)+
theme_cowplot()+
ylab("Simpson's index (inv.)")+
xlab("Infusion")+
scale_x_discrete(labels=c("OL","SG","PW"))+
scale_fill_manual(values=c("#B5D69F","#7DBA54","#C488A2","#992559","#B6F6EE","#60eede"),name="Infusion, temp.",labels=c("OL_cool","OL_warm","SG_cool","SG_warm","PW_cool","PW_warm"))+
ggtitle("Mesocosm water")
3.3 Statz
3.3.1 Richness
3.3.1.1 Water richness stats
df.div.mw <- subset(df.div.exp,type=="Meso. water")
df.div.mw$day <- as.factor(df.div.mw$day)
str(df.div.mw)## 'data.frame': 211 obs. of 30 variables:
## $ sample_name : chr "M_D13_OL1_1_S216_L001" "M_D13_OL1_2_S217_L001" "M_D13_OL1_3_S218_L001" "M_D13_OL1_4_S219_L001" ...
## $ Observed : num 25 28 22 28 31 31 28 29 27 21 ...
## $ Shannon : num 1.63 1.53 1.58 1.53 1.75 ...
## $ InvSimpson : num 2.67 2.4 2.44 2.32 2.93 ...
## $ orgname : chr "M_D13_OL1.1" "M_D13_OL1.2" "M_D13_OL1.3" "M_D13_OL1.4" ...
## $ mesocosm : chr "OL1.1" "OL1.2" "OL1.3" "OL1.4" ...
## $ type : Factor w/ 2 levels "Mosquitoes","Meso. water": 2 2 2 2 2 2 2 2 2 2 ...
## $ stage : chr "Water" "Water" "Water" "Water" ...
## $ sex : chr "" "" "" "" ...
## $ date.collected : chr "8/27/19" "8/27/19" "8/27/19" "8/27/19" ...
## $ day : Factor w/ 3 levels "4","12","20": 2 2 2 2 2 2 2 2 2 2 ...
## $ abdomen.length.mm: num NA NA NA NA NA NA NA NA NA NA ...
## $ wing.length : num NA NA NA NA NA NA NA NA NA NA ...
## $ hindleg.length : num NA NA NA NA NA NA NA NA NA NA ...
## $ special : chr "" "" "" "" ...
## $ dispersal : chr "N" "N" "D" "D" ...
## $ temperature : chr "C" "C" "C" "C" ...
## $ infusion : chr "OL" "OL" "OL" "OL" ...
## $ raw_miseq_reads : int 30391 39927 18829 58206 54860 39963 37857 45572 36841 27161 ...
## $ X16scq : num NA NA NA NA NA NA NA NA NA NA ...
## $ Wolb_ct1 : num NA NA NA NA NA NA NA NA NA NA ...
## $ Wolb_ct2 : num NA NA NA NA NA NA NA NA NA NA ...
## $ Act_ct : num NA NA NA NA NA NA NA NA NA NA ...
## $ notes : chr "" "" "" "" ...
## $ newday : chr "early" "early" "early" "early" ...
## $ lib_size : num 25797 34139 15861 49298 44559 ...
## $ is.neg : logi FALSE FALSE FALSE FALSE FALSE FALSE ...
## $ even : num 0.506 0.458 0.511 0.46 0.51 ...
## $ lib_size_trim : num 25797 34134 15861 49290 44524 ...
## $ inf.temp : chr "OLC" "OLC" "OLC" "OLC" ...#hist(df.div.mw$Observed)
#shapiro.test(log(df.div.mw$Observed))
#shapiro.test(df.div.mw$Observed)
wrich1<-glmmTMB(Observed~day+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),family="poisson",data=df.div.mw)
wrich2<-glmmTMB(Observed~day+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),family="compois",data=df.div.mw)
wrich3<-glmmTMB(Observed~day+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
AICtab(wrich1,wrich2,wrich3)## dAIC df
## wrich3 0.0 9
## wrich2 302.7 9
## wrich1 312.9 8##Gaussian best
##full model
wrich1<-glmmTMB(Observed~day+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
##no dispersal
wrich1nd<-glmmTMB(Observed~day+temperature+infusion+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mw)
##minus infusion
wrich1ni<-glmmTMB(Observed~day+temperature+dispersal+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mw)
##minus temperature
wrich1nt<-glmmTMB(Observed~day+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mw)
##minus time
wrich1nnd<-glmmTMB(Observed~temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mw)
anova(wrich1,wrich1nd) #0.3931## Data: df.div.mw
## Models:
## wrich1nd: Observed ~ day + temperature + infusion + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1nd: (1 | mesocosm), zi=~0, disp=~1
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1: (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1nd 8 1035.4 1062.2 -509.69 1019.4
## wrich1 9 1036.7 1066.8 -509.33 1018.7 0.7293 1 0.3931## Data: df.div.mw
## Models:
## wrich1ni: Observed ~ day + temperature + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1ni: (1 | mesocosm), zi=~0, disp=~1
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1: (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1ni 7 1118.0 1141.5 -552.02 1104.0
## wrich1 9 1036.7 1066.8 -509.33 1018.7 85.373 2 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: df.div.mw
## Models:
## wrich1nt: Observed ~ day + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1nt: (1 | mesocosm), zi=~0, disp=~1
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1: (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1nt 8 1035.3 1062.2 -509.67 1019.3
## wrich1 9 1036.7 1066.8 -509.33 1018.7 0.6732 1 0.4119## Data: df.div.mw
## Models:
## wrich1nnd: Observed ~ temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1nnd: (1 | mesocosm), zi=~0, disp=~1
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1: (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1nnd 7 1045.2 1068.7 -515.61 1031.2
## wrich1 9 1036.7 1066.8 -509.33 1018.7 12.559 2 0.001874 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Observed
## Chisq Df Pr(>Chisq)
## day 12.9404 2 0.001549 **
## temperature 0.6742 1 0.411578
## infusion 151.4172 2 < 2.2e-16 ***
## dispersal 0.7306 1 0.392704
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Response: Observed
# Chisq Df Pr(>Chisq)
# day 12.9404 2 0.001549 **
# temperature 0.6742 1 0.411577
# infusion 151.4172 2 < 2.2e-16 ***
# dispersal 0.7306 1 0.392704
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#AICtab(wrich1,wrich1nd,wrich1ni,wrich1nt,wrich1nnd) #no temp best... almost tied with no dispersal, then full
##interactions?
wrich1.int <- glmmTMB(Observed~day*infusion+temperature+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
wrich1.int1 <- glmmTMB(Observed~day+infusion*temperature+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
wrich1.int2 <- glmmTMB(Observed~day+infusion+temperature*dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
wrich1.int3 <- glmmTMB(Observed~day*temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
wrich1.int4 <- glmmTMB(Observed~day*dispersal+temperature+infusion+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
wrich1.int5 <- glmmTMB(Observed~day+dispersal*infusion+temperature+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
anova(wrich1,wrich1.int) #2.86e-10***## Data: df.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1: (1 | mesocosm), zi=~0, disp=~1
## wrich1.int: Observed ~ day * infusion + temperature + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1.int: (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1 9 1036.66 1066.8 -509.33 1018.66
## wrich1.int 13 994.79 1038.4 -484.40 968.79 49.863 4 3.857e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: df.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1: (1 | mesocosm), zi=~0, disp=~1
## wrich1.int1: Observed ~ day + infusion * temperature + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1.int1: (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1 9 1036.7 1066.8 -509.33 1018.7
## wrich1.int1 11 1037.2 1074.1 -507.62 1015.2 3.4154 2 0.1813## Data: df.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1: (1 | mesocosm), zi=~0, disp=~1
## wrich1.int2: Observed ~ day + infusion + temperature * dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1.int2: (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1 9 1036.7 1066.8 -509.33 1018.7
## wrich1.int2 10 1037.1 1070.7 -508.57 1017.1 1.5253 1 0.2168## Data: df.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1: (1 | mesocosm), zi=~0, disp=~1
## wrich1.int3: Observed ~ day * temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1.int3: (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1 9 1036.7 1066.8 -509.33 1018.7
## wrich1.int3 11 1039.0 1075.9 -508.52 1017.0 1.6226 2 0.4443## Data: df.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1: (1 | mesocosm), zi=~0, disp=~1
## wrich1.int4: Observed ~ day * dispersal + temperature + infusion + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1.int4: (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1 9 1036.7 1066.8 -509.33 1018.7
## wrich1.int4 11 1039.3 1076.2 -508.66 1017.3 1.3341 2 0.5132## Data: df.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1: (1 | mesocosm), zi=~0, disp=~1
## wrich1.int5: Observed ~ day + dispersal * infusion + temperature + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1.int5: (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1 9 1036.7 1066.8 -509.33 1018.7
## wrich1.int5 11 1037.8 1074.7 -507.90 1015.8 2.8599 2 0.2393## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Observed
## Chisq Df Pr(>Chisq)
## day 17.4815 2 0.0001599 ***
## infusion 167.8498 2 < 2.2e-16 ***
## temperature 0.7595 1 0.3834697
## dispersal 0.7626 1 0.3825074
## day:infusion 59.0003 4 4.705e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Analysis of Deviance Table (Type II Wald chisquare tests)
#
# Response: Observed
# Chisq Df Pr(>Chisq)
# day 17.4811 2 0.00016 ***
# infusion 167.8537 2 < 2.2e-16 ***
# temperature 0.7595 1 0.38348
# dispersal 0.7627 1 0.38248
# day:infusion 58.9948 4 4.717e-12 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#tapply(df.div.mw$log(Observed), df.div.mw$infusion, mean)
#tapply(df.div.mw$log(Observed), df.div.mw$day, mean)
wrich1.int.notem <- glmmTMB(Observed~day*infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
anova(wrich1.int,wrich1.int.notem) #0.38## Data: df.div.mw
## Models:
## wrich1.int.notem: Observed ~ day * infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1.int.notem: (1 | mesocosm), zi=~0, disp=~1
## wrich1.int: Observed ~ day * infusion + temperature + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1.int: (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1.int.notem 12 993.55 1033.8 -484.78 969.55
## wrich1.int 13 994.79 1038.4 -484.40 968.79 0.7567 1 0.3844wrich1.int.nodis <- glmmTMB(Observed~day*infusion+temperature+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
anova(wrich1.int,wrich1.int.nodis) #0.38## Data: df.div.mw
## Models:
## wrich1.int.nodis: Observed ~ day * infusion + temperature + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1.int.nodis: (1 | mesocosm), zi=~0, disp=~1
## wrich1.int: Observed ~ day * infusion + temperature + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1.int: (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1.int.nodis 12 993.55 1033.8 -484.78 969.55
## wrich1.int 13 994.79 1038.4 -484.40 968.79 0.7574 1 0.3841##
## Shapiro-Wilk normality test
##
## data: residuals(wrich1.int)
## W = 0.99201, p-value = 0.3043
##posthoc things
##just without dispersal because not plotting it
##also not plotting day but it should be in the model somewhere already
wrich1.int.em <- emmeans(wrich1.int,~infusion*temperature)
multcomp::cld(wrich1.int.em)## infusion temperature emmean SE df lower.CL upper.CL .group
## SW C 19.1 0.351 198 18.5 19.8 1
## SW H 19.5 0.353 198 18.8 20.1 1
## SG C 21.2 0.355 198 20.5 21.9 2
## SG H 21.5 0.357 198 20.8 22.2 2
## OL C 24.7 0.351 198 24.0 25.4 3
## OL H 25.0 0.353 198 24.3 25.7 3
##
## Results are averaged over the levels of: day, dispersal
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 6 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.3.3.1.2 Water richness stats - rarefied
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 102 taxa and 406 samples ]
## sample_data() Sample Data: [ 406 samples by 24 sample variables ]
## tax_table() Taxonomy Table: [ 102 taxa by 11 taxonomic ranks ]samdf.rare <- data.frame(ps.rare@sam_data)
df.rare <- data.frame(estimate_richness(ps.rare, split=TRUE, measures=c("Shannon","InvSimpson","Observed")))
df.rare$sample_name <- rownames(df.rare)
df.rare.div <- merge(df.rare,samdf.rare,by="sample_name") #add sample data
df.rare.div.exp <- subset(df.rare.div,type=="Microbial Water"|type=="A.albopictus")Same as above but without offset for library size
df.rare.div.mw <- subset(df.rare.div.exp,type=="Microbial Water")
df.rare.div.mw$day <- as.factor(df.rare.div.mw$day)
str(df.rare.div.mw)## 'data.frame': 211 obs. of 27 variables:
## $ sample_name : chr "M_D13_OL1_1_S216_L001" "M_D13_OL1_2_S217_L001" "M_D13_OL1_3_S218_L001" "M_D13_OL1_4_S219_L001" ...
## $ Observed : num 25 27 22 28 31 31 28 29 27 21 ...
## $ Shannon : num 1.64 1.51 1.56 1.5 1.77 ...
## $ InvSimpson : num 2.69 2.36 2.41 2.27 2.98 ...
## $ orgname : chr "M_D13_OL1.1" "M_D13_OL1.2" "M_D13_OL1.3" "M_D13_OL1.4" ...
## $ mesocosm : chr "OL1.1" "OL1.2" "OL1.3" "OL1.4" ...
## $ type : chr "Microbial Water" "Microbial Water" "Microbial Water" "Microbial Water" ...
## $ stage : chr "Water" "Water" "Water" "Water" ...
## $ sex : chr "" "" "" "" ...
## $ date.collected : chr "8/27/19" "8/27/19" "8/27/19" "8/27/19" ...
## $ day : Factor w/ 3 levels "4","12","20": 2 2 2 2 2 2 2 2 2 2 ...
## $ abdomen.length.mm: num NA NA NA NA NA NA NA NA NA NA ...
## $ wing.length : num NA NA NA NA NA NA NA NA NA NA ...
## $ hindleg.length : num NA NA NA NA NA NA NA NA NA NA ...
## $ special : chr "" "" "" "" ...
## $ dispersal : chr "N" "N" "D" "D" ...
## $ temperature : chr "C" "C" "C" "C" ...
## $ infusion : chr "OL" "OL" "OL" "OL" ...
## $ raw_miseq_reads : int 30391 39927 18829 58206 54860 39963 37857 45572 36841 27161 ...
## $ X16scq : num NA NA NA NA NA NA NA NA NA NA ...
## $ Wolb_ct1 : num NA NA NA NA NA NA NA NA NA NA ...
## $ Wolb_ct2 : num NA NA NA NA NA NA NA NA NA NA ...
## $ Act_ct : num NA NA NA NA NA NA NA NA NA NA ...
## $ notes : chr "" "" "" "" ...
## $ newday : chr "early" "early" "early" "early" ...
## $ lib_size : num 25797 34139 15861 49298 44559 ...
## $ is.neg : logi FALSE FALSE FALSE FALSE FALSE FALSE ...#hist(df.rare.div.mw$Observed)
#shapiro.test(log(df.rare.div.mw$Observed))
#shapiro.test(df.rare.div.mw$Observed)
#wrich1<-glmmTMB(Observed~day+temperature+infusion+dispersal+(1|mesocosm),family="poisson",data=df.rare.div.mw)
#wrich2<-glmmTMB(Observed~day+temperature+infusion+dispersal+(1|mesocosm),family="compois",data=df.rare.div.mw)
wrich3<-glmmTMB(Observed~day+temperature+infusion+dispersal+(1|mesocosm),data=df.rare.div.mw)
#AICtab(wrich1,wrich2,wrich3)
##full model
wrich1<-glmmTMB(Observed~day+temperature+infusion+dispersal+(1|mesocosm),data=df.rare.div.mw)
##no dispersal
wrich1nd<-glmmTMB(Observed~day+temperature+infusion+(1|mesocosm), data=df.rare.div.mw)
##minus infusion
wrich1ni<-glmmTMB(Observed~day+temperature+dispersal+(1|mesocosm), data=df.rare.div.mw)
##minus temperature
wrich1nt<-glmmTMB(Observed~day+infusion+dispersal+(1|mesocosm), data=df.rare.div.mw)
##minus time
wrich1nnd<-glmmTMB(Observed~temperature+infusion+dispersal+(1|mesocosm), data=df.rare.div.mw)
anova(wrich1,wrich1nd) #0.466## Data: df.rare.div.mw
## Models:
## wrich1nd: Observed ~ day + temperature + infusion + (1 | mesocosm), zi=~0, disp=~1
## wrich1: Observed ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1nd 8 1066.8 1093.7 -525.43 1050.8
## wrich1 9 1068.3 1098.5 -525.16 1050.3 0.5315 1 0.466## Data: df.rare.div.mw
## Models:
## wrich1ni: Observed ~ day + temperature + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wrich1: Observed ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1ni 7 1145.2 1168.7 -565.62 1131.2
## wrich1 9 1068.3 1098.5 -525.16 1050.3 80.917 2 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: df.rare.div.mw
## Models:
## wrich1nt: Observed ~ day + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wrich1: Observed ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1nt 8 1066.8 1093.6 -525.41 1050.8
## wrich1 9 1068.3 1098.5 -525.16 1050.3 0.5014 1 0.4789## Data: df.rare.div.mw
## Models:
## wrich1nnd: Observed ~ temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wrich1: Observed ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1nnd 7 1074.3 1097.8 -530.17 1060.3
## wrich1 9 1068.3 1098.5 -525.16 1050.3 10.013 2 0.006693 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Observed
## Chisq Df Pr(>Chisq)
## day 10.2547 2 0.005932 **
## temperature 0.5020 1 0.478603
## infusion 133.7700 2 < 2.2e-16 ***
## dispersal 0.5322 1 0.465696
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Analysis of Deviance Table (Type II Wald chisquare tests)
#
# Response: Observed
# Chisq Df Pr(>Chisq)
# day 10.2547 2 0.005932 **
# temperature 0.5020 1 0.478604
# infusion 133.7700 2 < 2.2e-16 ***
# dispersal 0.5322 1 0.465695
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#AICtab(wrich1,wrich1nd,wrich1ni,wrich1nt,wrich1nnd) #no temp best... almost tied with no dispersal, then full
##interactions?
wrich1.int <- glmmTMB(Observed~day*infusion+temperature+dispersal+(1|mesocosm),data=df.rare.div.mw)
wrich1.int1 <- glmmTMB(Observed~day+infusion*temperature+dispersal+(1|mesocosm),data=df.rare.div.mw)
wrich1.int2 <- glmmTMB(Observed~day+infusion+temperature*dispersal+(1|mesocosm),data=df.rare.div.mw)
wrich1.int3 <- glmmTMB(Observed~day*temperature+infusion+dispersal+(1|mesocosm),data=df.rare.div.mw)
wrich1.int4 <- glmmTMB(Observed~day*dispersal+temperature+infusion+(1|mesocosm),data=df.rare.div.mw)
wrich1.int5 <- glmmTMB(Observed~day+dispersal*infusion+temperature+(1|mesocosm),data=df.rare.div.mw)
anova(wrich1,wrich1.int) #4.35e-11***## Data: df.rare.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wrich1.int: Observed ~ day * infusion + temperature + dispersal + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1 9 1068.3 1098.5 -525.16 1050.32
## wrich1.int 13 1021.9 1065.5 -497.96 995.93 54.391 4 4.358e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: df.rare.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wrich1.int1: Observed ~ day + infusion * temperature + dispersal + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1 9 1068.3 1098.5 -525.16 1050.3
## wrich1.int1 11 1068.7 1105.5 -523.33 1046.7 3.6606 2 0.1604## Data: df.rare.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wrich1.int2: Observed ~ day + infusion + temperature * dispersal + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1 9 1068.3 1098.5 -525.16 1050.3
## wrich1.int2 10 1068.5 1102.0 -524.24 1048.5 1.836 1 0.1754## Data: df.rare.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wrich1.int3: Observed ~ day * temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1 9 1068.3 1098.5 -525.16 1050.3
## wrich1.int3 11 1071.0 1107.8 -524.48 1049.0 1.3561 2 0.5076## Data: df.rare.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wrich1.int4: Observed ~ day * dispersal + temperature + infusion + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1 9 1068.3 1098.5 -525.16 1050.3
## wrich1.int4 11 1071.5 1108.3 -524.74 1049.5 0.8411 2 0.6567## Data: df.rare.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wrich1.int5: Observed ~ day + dispersal * infusion + temperature + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1 9 1068.3 1098.5 -525.16 1050.3
## wrich1.int5 11 1070.2 1107.1 -524.10 1048.2 2.1147 2 0.3474## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Observed
## Chisq Df Pr(>Chisq)
## day 13.9846 2 0.0009189 ***
## infusion 154.5864 2 < 2.2e-16 ***
## temperature 0.5795 1 0.4464952
## dispersal 0.5809 1 0.4459485
## day:infusion 64.8561 4 2.759e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Analysis of Deviance Table (Type II Wald chisquare tests)
#
# Response: Observed
# Chisq Df Pr(>Chisq)
# day 13.9847 2 0.0009189 ***
# infusion 154.5860 2 < 2.2e-16 ***
# temperature 0.5795 1 0.4464988
# dispersal 0.5809 1 0.4459484
# day:infusion 64.8562 4 2.759e-13 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
wrich1.int.notem <- glmmTMB(Observed~day*infusion+dispersal+(1|mesocosm),data=df.rare.div.mw)
anova(wrich1.int,wrich1.int.notem) #0.45## Data: df.rare.div.mw
## Models:
## wrich1.int.notem: Observed ~ day * infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wrich1.int: Observed ~ day * infusion + temperature + dispersal + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1.int.notem 12 1020.5 1060.7 -498.25 996.51
## wrich1.int 13 1021.9 1065.5 -497.96 995.93 0.5777 1 0.4472wrich1.int.nodis <- glmmTMB(Observed~day*infusion+temperature+(1|mesocosm),data=df.rare.div.mw)
anova(wrich1.int,wrich1.int.nodis) #0.45## Data: df.rare.div.mw
## Models:
## wrich1.int.nodis: Observed ~ day * infusion + temperature + (1 | mesocosm), zi=~0, disp=~1
## wrich1.int: Observed ~ day * infusion + temperature + dispersal + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1.int.nodis 12 1020.5 1060.7 -498.25 996.51
## wrich1.int 13 1021.9 1065.5 -497.96 995.93 0.5778 1 0.4472#tapply(df.rare.div.mw$log(Observed), df.rare.div.mw$infusion, mean)
#tapply(df.rare.div.mw$log(Observed), df.rare.div.mw$day, mean)
##model checking
shapiro.test(residuals(wrich1.int)) #awesome##
## Shapiro-Wilk normality test
##
## data: residuals(wrich1.int)
## W = 0.99017, p-value = 0.1609
##posthoc things
##just without dispersal because not plotting it
##also not plotting day but it should be in the model somewhere already
wrich1.int.em <- emmeans(wrich1.int,~infusion*temperature)
multcomp::cld(wrich1.int.em)## infusion temperature emmean SE df lower.CL upper.CL .group
## SW C 19.1 0.371 198 18.3 19.8 1
## SW H 19.4 0.373 198 18.6 20.1 12
## SG C 20.9 0.374 198 20.1 21.6 23
## SG H 21.1 0.377 198 20.4 21.9 3
## OL C 24.6 0.371 198 23.9 25.4 4
## OL H 24.9 0.373 198 24.2 25.7 4
##
## Results are averaged over the levels of: day, dispersal
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 6 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.3.3.1.3 Mosquito richness stats
##
## Shapiro-Wilk normality test
##
## data: df.div.mq$Observed
## W = 0.89464, p-value = 1.697e-10##different stat families
mrich1 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
mrich2 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),family="compois",data=df.div.mq)
mrich3 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),family="poisson",data=df.div.mq)
AICtab(mrich1,mrich2,mrich3) #gaussian again## dAIC df
## mrich1 0.0 10
## mrich2 47.1 10
## mrich3 48.3 9##full model
mrich1 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
##no dispersal
mrich1nd <- glmmTMB(Observed~newday+temperature+infusion+sex+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
##minus infusion
mrich1ni <- glmmTMB(Observed~newday+temperature+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
##minus temperature
mrich1nt <- glmmTMB(Observed~newday+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
##minus time
mrich1nnd <- glmmTMB(Observed~temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mq)
##minus sex
mrich1ns <- glmmTMB(Observed~newday+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
anova(mrich1,mrich1nd) #0.24## Data: df.div.mq
## Models:
## mrich1nd: Observed ~ newday + temperature + infusion + sex + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## mrich1nd: (1 | mesocosm), zi=~0, disp=~1
## mrich1: Observed ~ newday + temperature + infusion + dispersal + sex + , zi=~0, disp=~1
## mrich1: offset(log(lib_size_trim)) + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mrich1nd 9 917.87 947.32 -449.93 899.87
## mrich1 10 918.46 951.19 -449.23 898.46 1.4056 1 0.2358## Data: df.div.mq
## Models:
## mrich1ni: Observed ~ newday + temperature + dispersal + sex + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## mrich1ni: (1 | mesocosm), zi=~0, disp=~1
## mrich1: Observed ~ newday + temperature + infusion + dispersal + sex + , zi=~0, disp=~1
## mrich1: offset(log(lib_size_trim)) + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mrich1ni 8 914.95 941.14 -449.48 898.95
## mrich1 10 918.46 951.19 -449.23 898.46 0.4942 2 0.7811## Data: df.div.mq
## Models:
## mrich1nt: Observed ~ newday + infusion + dispersal + sex + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## mrich1nt: (1 | mesocosm), zi=~0, disp=~1
## mrich1: Observed ~ newday + temperature + infusion + dispersal + sex + , zi=~0, disp=~1
## mrich1: offset(log(lib_size_trim)) + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mrich1nt 9 916.68 946.14 -449.34 898.68
## mrich1 10 918.46 951.19 -449.23 898.46 0.2243 1 0.6358## Data: df.div.mq
## Models:
## mrich1nnd: Observed ~ temperature + infusion + dispersal + sex + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## mrich1nnd: (1 | mesocosm), zi=~0, disp=~1
## mrich1: Observed ~ newday + temperature + infusion + dispersal + sex + , zi=~0, disp=~1
## mrich1: offset(log(lib_size_trim)) + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mrich1nnd 8 916.27 942.46 -450.14 900.27
## mrich1 10 918.46 951.19 -449.23 898.46 1.814 2 0.4037## Data: df.div.mq
## Models:
## mrich1ns: Observed ~ newday + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## mrich1ns: (1 | mesocosm), zi=~0, disp=~1
## mrich1: Observed ~ newday + temperature + infusion + dispersal + sex + , zi=~0, disp=~1
## mrich1: offset(log(lib_size_trim)) + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mrich1ns 9 916.86 946.32 -449.43 898.86
## mrich1 10 918.46 951.19 -449.23 898.46 0.3989 1 0.5277## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Observed
## Chisq Df Pr(>Chisq)
## newday 1.8273 2 0.4011
## temperature 0.2244 1 0.6357
## infusion 0.4926 2 0.7817
## dispersal 1.4646 1 0.2262
## sex 0.3992 1 0.5275##
## Shapiro-Wilk normality test
##
## data: residuals(mrich1)
## W = 0.90866, p-value = 1.328e-09
## infusion temperature emmean SE df lower.CL upper.CL .group
## OL H 8.34 0.417 185 7.51 9.16 1
## SG H 8.49 0.447 185 7.61 9.37 1
## OL C 8.56 0.340 185 7.89 9.23 1
## SG C 8.72 0.490 185 7.75 9.68 1
## SW H 8.77 0.479 185 7.82 9.71 1
## SW C 8.99 0.564 185 7.88 10.10 1
##
## Results are averaged over the levels of: newday, dispersal, sex
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 6 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.3.3.1.4 Mosquito richness stats - rarefied
Same as above but without offset for library size
df.rare.div.mq <- subset(df.rare.div.exp,type=="A.albopictus")
#hist(df.rare.div.mq$Observed)
#shapiro.test(log(df.rare.div.mq$Observed))
#shapiro.test(df.rare.div.mq$Observed)
#mrich1<-glmmTMB(Observed~day+temperature+infusion+dispersal+(1|mesocosm),family="poisson",data=df.rare.div.mq)
#mrich2<-glmmTMB(Observed~day+temperature+infusion+dispersal+(1|mesocosm),family="compois",data=df.rare.div.mq)
#mrich3<-glmmTMB(Observed~day+temperature+infusion+dispersal+(1|mesocosm),data=df.rare.div.mq)
#AICtab(mrich1,mrich2,mrich3)
##full model
mrich1<-glmmTMB(Observed~newday+temperature+infusion+dispersal+(1|mesocosm),data=df.rare.div.mq)
##no dispersal
mrich1nd<-glmmTMB(Observed~newday+temperature+infusion+(1|mesocosm), data=df.rare.div.mq)
##minus infusion
mrich1ni<-glmmTMB(Observed~newday+temperature+dispersal+(1|mesocosm), data=df.rare.div.mq)
##minus temperature
mrich1nt<-glmmTMB(Observed~newday+infusion+dispersal+(1|mesocosm), data=df.rare.div.mq)
##minus time
mrich1nnd<-glmmTMB(Observed~temperature+infusion+dispersal+(1|mesocosm), data=df.rare.div.mq)
anova(mrich1,mrich1nd) #0.18## Data: df.rare.div.mq
## Models:
## mrich1nd: Observed ~ newday + temperature + infusion + (1 | mesocosm), zi=~0, disp=~1
## mrich1: Observed ~ newday + temperature + infusion + dispersal + (1 | , zi=~0, disp=~1
## mrich1: mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mrich1nd 8 914.40 940.58 -449.2 898.40
## mrich1 9 914.61 944.06 -448.3 896.61 1.7942 1 0.1804## Data: df.rare.div.mq
## Models:
## mrich1ni: Observed ~ newday + temperature + dispersal + (1 | mesocosm), zi=~0, disp=~1
## mrich1: Observed ~ newday + temperature + infusion + dispersal + (1 | , zi=~0, disp=~1
## mrich1: mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mrich1ni 7 910.75 933.66 -448.37 896.75
## mrich1 9 914.61 944.06 -448.30 896.61 0.144 2 0.9305## Data: df.rare.div.mq
## Models:
## mrich1nt: Observed ~ newday + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## mrich1: Observed ~ newday + temperature + infusion + dispersal + (1 | , zi=~0, disp=~1
## mrich1: mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mrich1nt 8 912.67 938.85 -448.33 896.67
## mrich1 9 914.61 944.06 -448.30 896.61 0.0605 1 0.8056## Data: df.rare.div.mq
## Models:
## mrich1nnd: Observed ~ temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## mrich1: Observed ~ newday + temperature + infusion + dispersal + (1 | , zi=~0, disp=~1
## mrich1: mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mrich1nnd 7 913.23 936.14 -449.61 899.23
## mrich1 9 914.61 944.06 -448.30 896.61 2.6238 2 0.2693## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Observed
## Chisq Df Pr(>Chisq)
## newday 2.7059 2 0.2585
## temperature 0.0606 1 0.8056
## infusion 0.1441 2 0.9305
## dispersal 1.8600 1 0.1726# Analysis of Deviance Table (Type II Wald chisquare tests)
#
# Response: Observed
# Chisq Df Pr(>Chisq)
# newday 2.7059 2 0.2585
# temperature 0.0606 1 0.8056
# infusion 0.1441 2 0.9305
# dispersal 1.8600 1 0.1726
#AICtab(mrich1,mrich1nd,mrich1ni,mrich1nt,mrich1nnd) #no temp best... almost tied with no dispersal, then full
##interactions?
mrich1.int <- glmmTMB(Observed~newday*infusion+temperature+dispersal+(1|mesocosm),data=df.rare.div.mq)## dropping columns from rank-deficient conditional model: newdaymid:infusionSWmrich1.int1 <- glmmTMB(Observed~newday+infusion*temperature+dispersal+(1|mesocosm),data=df.rare.div.mq)
mrich1.int2 <- glmmTMB(Observed~newday+infusion+temperature*dispersal+(1|mesocosm),data=df.rare.div.mq)
mrich1.int3 <- glmmTMB(Observed~newday*temperature+infusion+dispersal+(1|mesocosm),data=df.rare.div.mq)## dropping columns from rank-deficient conditional model: newdaymid:temperatureHmrich1.int4 <- glmmTMB(Observed~newday*dispersal+temperature+infusion+(1|mesocosm),data=df.rare.div.mq)
mrich1.int5 <- glmmTMB(Observed~newday+dispersal*infusion+temperature+(1|mesocosm),data=df.rare.div.mq)
anova(mrich1,mrich1.int) #ns## Data: df.rare.div.mq
## Models:
## mrich1: Observed ~ newday + temperature + infusion + dispersal + (1 | , zi=~0, disp=~1
## mrich1: mesocosm), zi=~0, disp=~1
## mrich1.int: Observed ~ newday * infusion + temperature + dispersal + (1 | , zi=~0, disp=~1
## mrich1.int: mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mrich1 9 914.61 944.06 -448.30 896.61
## mrich1.int 12 918.16 957.43 -447.08 894.16 2.4499 3 0.4844## Data: df.rare.div.mq
## Models:
## mrich1: Observed ~ newday + temperature + infusion + dispersal + (1 | , zi=~0, disp=~1
## mrich1: mesocosm), zi=~0, disp=~1
## mrich1.int1: Observed ~ newday + infusion * temperature + dispersal + (1 | , zi=~0, disp=~1
## mrich1.int1: mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mrich1 9 914.61 944.06 -448.30 896.61
## mrich1.int1 11 918.12 954.12 -448.06 896.12 0.4879 2 0.7835## Data: df.rare.div.mq
## Models:
## mrich1: Observed ~ newday + temperature + infusion + dispersal + (1 | , zi=~0, disp=~1
## mrich1: mesocosm), zi=~0, disp=~1
## mrich1.int2: Observed ~ newday + infusion + temperature * dispersal + (1 | , zi=~0, disp=~1
## mrich1.int2: mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mrich1 9 914.61 944.06 -448.30 896.61
## mrich1.int2 10 915.93 948.66 -447.97 895.93 0.6735 1 0.4118## Data: df.rare.div.mq
## Models:
## mrich1: Observed ~ newday + temperature + infusion + dispersal + (1 | , zi=~0, disp=~1
## mrich1: mesocosm), zi=~0, disp=~1
## mrich1.int3: Observed ~ newday * temperature + infusion + dispersal + (1 | , zi=~0, disp=~1
## mrich1.int3: mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mrich1 9 914.61 944.06 -448.30 896.61
## mrich1.int3 10 916.54 949.27 -448.27 896.54 0.0623 1 0.8029## Data: df.rare.div.mq
## Models:
## mrich1: Observed ~ newday + temperature + infusion + dispersal + (1 | , zi=~0, disp=~1
## mrich1: mesocosm), zi=~0, disp=~1
## mrich1.int4: Observed ~ newday * dispersal + temperature + infusion + (1 | , zi=~0, disp=~1
## mrich1.int4: mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mrich1 9 914.61 944.06 -448.30 896.61
## mrich1.int4 11 917.74 953.74 -447.87 895.74 0.8642 2 0.6491## Data: df.rare.div.mq
## Models:
## mrich1: Observed ~ newday + temperature + infusion + dispersal + (1 | , zi=~0, disp=~1
## mrich1: mesocosm), zi=~0, disp=~1
## mrich1.int5: Observed ~ newday + dispersal * infusion + temperature + (1 | , zi=~0, disp=~1
## mrich1.int5: mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mrich1 9 914.61 944.06 -448.30 896.61
## mrich1.int5 11 916.84 952.84 -447.42 894.84 1.7644 2 0.4139## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Observed
## Chisq Df Pr(>Chisq)
## newday 2.7402 2 0.2541
## infusion 0.1459 2 0.9297
## temperature 0.0121 1 0.9125
## dispersal 2.4306 1 0.1190
## newday:infusion 2.4653 3 0.4816# Response: Observed
# Chisq Df Pr(>Chisq)
# newday 2.7401 2 0.2541
# infusion 0.1459 2 0.9296
# temperature 0.0121 1 0.9125
# dispersal 2.4306 1 0.1190
# newday:infusion 2.4653 3 0.4816
#tapply(df.rare.div.mq$log(Observed), df.rare.div.mq$infusion, mean)
#tapply(df.rare.div.mq$log(Observed), df.rare.div.mq$day, mean)
##model checking
shapiro.test(residuals(mrich1))##
## Shapiro-Wilk normality test
##
## data: residuals(mrich1)
## W = 0.90715, p-value = 1.054e-09
##posthoc things
##just without dispersal because not plotting it
##also not plotting day but it should be in the model somewhere already
mrich1.em <- emmeans(mrich1,~infusion+temperature)
multcomp::cld(mrich1.em)## infusion temperature emmean SE df lower.CL upper.CL .group
## OL H 7.91 0.410 186 7.10 8.72 1
## SG H 8.02 0.435 186 7.16 8.88 1
## OL C 8.03 0.328 186 7.38 8.67 1
## SW H 8.13 0.463 186 7.22 9.04 1
## SG C 8.14 0.478 186 7.19 9.08 1
## SW C 8.24 0.552 186 7.16 9.33 1
##
## Results are averaged over the levels of: newday, dispersal
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 6 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.3.3.2 Simpson’s
3.3.2.1 Water Simpson’s stats
##
## Shapiro-Wilk normality test
##
## data: df.div.mw$InvSimpson
## W = 0.96756, p-value = 9.003e-05
##
## Shapiro-Wilk normality test
##
## data: log(df.div.mw$InvSimpson)
## W = 0.95673, p-value = 5.166e-06## Best Normalizing transformation with 211 Observations
## Estimated Normality Statistics (Pearson P / df, lower => more normal):
## - arcsinh(x): 1.5509
## - Box-Cox: 1.4324
## - Center+scale: 1.4561
## - Double Reversed Log_b(x+a): 1.6847
## - Exp(x): 3.3273
## - Log_b(x+a): 1.6108
## - orderNorm (ORQ): 1.1559
## - sqrt(x + a): 1.4286
## - Yeo-Johnson: 1.4424
## Estimation method: Out-of-sample via CV with 10 folds and 5 repeats
##
## Based off these, bestNormalize chose:
## orderNorm Transformation with 211 nonmissing obs and no ties
## - Original quantiles:
## 0% 25% 50% 75% 100%
## 1.359 1.947 2.669 3.225 5.130##
## Shapiro-Wilk normality test
##
## data: df.div.mw$simp.norm
## W = 0.99978, p-value = 1##just replaced InvSimpson with simp.norm
##full model
wsimp1 <- glmmTMB(simp.norm~day+temperature+infusion+dispersal+(1|mesocosm),data=df.div.mw)
##no dispersal
wsimp1nd <- glmmTMB(simp.norm~day+temperature+infusion+(1|mesocosm),data=df.div.mw)
##minus infusion
wsimp1ni <- glmmTMB(simp.norm~day+temperature+dispersal+(1|mesocosm),data=df.div.mw)
##minus temperature
wsimp1nt <- glmmTMB(simp.norm~day+infusion+dispersal+(1|mesocosm),data=df.div.mw)
##minus time
wsimp1nnd <- glmmTMB(simp.norm~temperature+infusion+dispersal+(1|mesocosm),data=df.div.mw)
anova(wsimp1,wsimp1nd) #0.86## Data: df.div.mw
## Models:
## wsimp1nd: simp.norm ~ day + temperature + infusion + (1 | mesocosm), zi=~0, disp=~1
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wsimp1nd 8 388.09 414.90 -186.04 372.09
## wsimp1 9 390.08 420.24 -186.04 372.08 0.0091 1 0.924## Data: df.div.mw
## Models:
## wsimp1ni: simp.norm ~ day + temperature + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wsimp1ni 7 489.84 513.30 -237.92 475.84
## wsimp1 9 390.08 420.24 -186.04 372.08 103.76 2 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: df.div.mw
## Models:
## wsimp1nt: simp.norm ~ day + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wsimp1nt 8 420.66 447.48 -202.33 404.66
## wsimp1 9 390.08 420.24 -186.04 372.08 32.585 1 1.141e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: df.div.mw
## Models:
## wsimp1nnd: simp.norm ~ temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wsimp1nnd 7 428.81 452.27 -207.40 414.81
## wsimp1 9 390.08 420.24 -186.04 372.08 42.732 2 5.259e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: simp.norm
## Chisq Df Pr(>Chisq)
## day 49.7901 2 1.542e-11 ***
## temperature 41.3001 1 1.306e-10 ***
## infusion 229.1056 2 < 2.2e-16 ***
## dispersal 0.0091 1 0.924
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Analysis of Deviance Table (Type II Wald chisquare tests)
#
# Response: simp.norm
# Chisq Df Pr(>Chisq)
# day 45.4454 2 1.354e-10 ***
# temperature 48.9643 1 2.607e-12 ***
# infusion 276.6117 2 < 2.2e-16 ***
# dispersal 0.0291 1 0.8644
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##interaction?
wsimp1.int <- glmmTMB(simp.norm~day*infusion+temperature+dispersal+(1|mesocosm),data=df.div.mw)
wsimp1.int1 <- glmmTMB(simp.norm~day+infusion*temperature+dispersal+(1|mesocosm),data=df.div.mw)
wsimp1.int2 <- glmmTMB(simp.norm~day+infusion+temperature*dispersal+(1|mesocosm),data=df.div.mw)
wsimp1.int3 <- glmmTMB(simp.norm~day*temperature+infusion+dispersal+(1|mesocosm),data=df.div.mw)
wsimp1.int4 <- glmmTMB(simp.norm~day*dispersal+temperature+infusion+(1|mesocosm),data=df.div.mw)
wsimp1.int5 <- glmmTMB(simp.norm~day+dispersal*infusion+temperature+(1|mesocosm),data=df.div.mw)
anova(wsimp1,wsimp1.int) #4.27e-09***## Data: df.div.mw
## Models:
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1.int: simp.norm ~ day * infusion + temperature + dispersal + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wsimp1 9 390.08 420.24 -186.04 372.08
## wsimp1.int 13 349.46 393.04 -161.73 323.46 48.615 4 7.026e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: df.div.mw
## Models:
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1.int1: simp.norm ~ day + infusion * temperature + dispersal + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wsimp1 9 390.08 420.24 -186.04 372.08
## wsimp1.int1 11 389.97 426.84 -183.98 367.97 4.1096 2 0.1281## Data: df.div.mw
## Models:
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1.int2: simp.norm ~ day + infusion + temperature * dispersal + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wsimp1 9 390.08 420.24 -186.04 372.08
## wsimp1.int2 10 391.39 424.90 -185.69 371.39 0.6898 1 0.4062## Data: df.div.mw
## Models:
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1.int3: simp.norm ~ day * temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wsimp1 9 390.08 420.24 -186.04 372.08
## wsimp1.int3 11 393.77 430.64 -185.89 371.77 0.3052 2 0.8585## Data: df.div.mw
## Models:
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1.int4: simp.norm ~ day * dispersal + temperature + infusion + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wsimp1 9 390.08 420.24 -186.04 372.08
## wsimp1.int4 11 392.09 428.96 -185.04 370.09 1.9908 2 0.3696## Data: df.div.mw
## Models:
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1.int5: simp.norm ~ day + dispersal * infusion + temperature + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wsimp1 9 390.08 420.24 -186.04 372.08
## wsimp1.int5 11 393.48 430.35 -185.74 371.48 0.5976 2 0.7417## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: simp.norm
## Chisq Df Pr(>Chisq)
## day 70.7469 2 4.340e-16 ***
## infusion 227.4983 2 < 2.2e-16 ***
## temperature 39.6143 1 3.094e-10 ***
## dispersal 0.0146 1 0.904
## day:infusion 58.5469 4 5.858e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Response: simp.norm
# Chisq Df Pr(>Chisq)
# day 62.7925 2 2.316e-14 ***
# infusion 274.8962 2 < 2.2e-16 ***
# temperature 47.2450 1 6.264e-12 ***
# dispersal 0.0397 1 0.842
# day:infusion 53.2218 4 7.658e-11 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
wsimp1.int.notem <- glmmTMB(simp.norm~day*infusion+dispersal+(1|mesocosm),data=df.div.mw)
anova(wsimp1.int,wsimp1.int.notem) #sig***## Data: df.div.mw
## Models:
## wsimp1.int.notem: simp.norm ~ day * infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1.int: simp.norm ~ day * infusion + temperature + dispersal + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wsimp1.int.notem 12 378.91 419.13 -177.46 354.91
## wsimp1.int 13 349.46 393.04 -161.73 323.46 31.451 1 2.046e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1wsimp1.int.nodis <- glmmTMB(simp.norm~day*infusion+temperature+(1|mesocosm),data=df.div.mw)
anova(wsimp1.int,wsimp1.int.nodis) #0.84## Data: df.div.mw
## Models:
## wsimp1.int.nodis: simp.norm ~ day * infusion + temperature + (1 | mesocosm), zi=~0, disp=~1
## wsimp1.int: simp.norm ~ day * infusion + temperature + dispersal + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wsimp1.int.nodis 12 347.48 387.70 -161.74 323.48
## wsimp1.int 13 349.46 393.04 -161.73 323.46 0.0146 1 0.904##
## Shapiro-Wilk normality test
##
## data: residuals(wsimp1.int)
## W = 0.93884, p-value = 9.519e-08
## infusion temperature emmean SE df lower.CL upper.CL .group
## SG C -1.2595 0.0939 198 -1.445 -1.0743 1
## SG H -0.6695 0.0943 198 -0.855 -0.4835 2
## OL C -0.0999 0.0934 198 -0.284 0.0842 3
## SW C 0.4404 0.0934 198 0.256 0.6245 4
## OL H 0.4901 0.0937 198 0.305 0.6750 4
## SW H 1.0304 0.0937 198 0.846 1.2153 5
##
## Results are averaged over the levels of: day, dispersal
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 6 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same. # infusion temperature emmean SE df lower.CL upper.CL .group
# SG C -1.300 0.0886 198 -1.474 -1.1249 1
# SG H -0.692 0.0890 198 -0.867 -0.5162 2
# OL C -0.105 0.0881 198 -0.278 0.0691 3
# SW C 0.463 0.0881 198 0.289 0.6365 4
# OL H 0.503 0.0884 198 0.329 0.6777 4
# SW H 1.071 0.0884 198 0.896 1.2452 53.3.2.2 Water Simpson’s stats - rarefied
## 'data.frame': 211 obs. of 27 variables:
## $ sample_name : chr "M_D13_OL1_1_S216_L001" "M_D13_OL1_2_S217_L001" "M_D13_OL1_3_S218_L001" "M_D13_OL1_4_S219_L001" ...
## $ Observed : num 25 27 22 28 31 31 28 29 27 21 ...
## $ Shannon : num 1.64 1.51 1.56 1.5 1.77 ...
## $ InvSimpson : num 2.69 2.36 2.41 2.27 2.98 ...
## $ orgname : chr "M_D13_OL1.1" "M_D13_OL1.2" "M_D13_OL1.3" "M_D13_OL1.4" ...
## $ mesocosm : chr "OL1.1" "OL1.2" "OL1.3" "OL1.4" ...
## $ type : chr "Microbial Water" "Microbial Water" "Microbial Water" "Microbial Water" ...
## $ stage : chr "Water" "Water" "Water" "Water" ...
## $ sex : chr "" "" "" "" ...
## $ date.collected : chr "8/27/19" "8/27/19" "8/27/19" "8/27/19" ...
## $ day : Factor w/ 3 levels "4","12","20": 2 2 2 2 2 2 2 2 2 2 ...
## $ abdomen.length.mm: num NA NA NA NA NA NA NA NA NA NA ...
## $ wing.length : num NA NA NA NA NA NA NA NA NA NA ...
## $ hindleg.length : num NA NA NA NA NA NA NA NA NA NA ...
## $ special : chr "" "" "" "" ...
## $ dispersal : chr "N" "N" "D" "D" ...
## $ temperature : chr "C" "C" "C" "C" ...
## $ infusion : chr "OL" "OL" "OL" "OL" ...
## $ raw_miseq_reads : int 30391 39927 18829 58206 54860 39963 37857 45572 36841 27161 ...
## $ X16scq : num NA NA NA NA NA NA NA NA NA NA ...
## $ Wolb_ct1 : num NA NA NA NA NA NA NA NA NA NA ...
## $ Wolb_ct2 : num NA NA NA NA NA NA NA NA NA NA ...
## $ Act_ct : num NA NA NA NA NA NA NA NA NA NA ...
## $ notes : chr "" "" "" "" ...
## $ newday : chr "early" "early" "early" "early" ...
## $ lib_size : num 25797 34139 15861 49298 44559 ...
## $ is.neg : logi FALSE FALSE FALSE FALSE FALSE FALSE ...## Best Normalizing transformation with 211 Observations
## Estimated Normality Statistics (Pearson P / df, lower => more normal):
## - arcsinh(x): 1.697
## - Box-Cox: 1.6111
## - Center+scale: 1.5965
## - Double Reversed Log_b(x+a): 1.5648
## - Exp(x): 3.1941
## - Log_b(x+a): 1.6639
## - orderNorm (ORQ): 1.3294
## - sqrt(x + a): 1.5573
## - Yeo-Johnson: 1.5903
## Estimation method: Out-of-sample via CV with 10 folds and 5 repeats
##
## Based off these, bestNormalize chose:
## orderNorm Transformation with 211 nonmissing obs and no ties
## - Original quantiles:
## 0% 25% 50% 75% 100%
## 1.367 1.957 2.684 3.244 5.127df.rare.div.mw$simp.norm <- bestNormalize(df.rare.div.mw$InvSimpson)$x.t
shapiro.test(df.rare.div.mw$simp.norm)##
## Shapiro-Wilk normality test
##
## data: df.rare.div.mw$simp.norm
## W = 0.99978, p-value = 1#hist(df.rare.div.mw$Observed)
#shapiro.test(log(df.rare.div.mw$Observed))
#shapiro.test(df.rare.div.mw$Observed)
##full model
wsimp1<-glmmTMB(simp.norm~day+temperature+infusion+dispersal+(1|mesocosm),data=df.rare.div.mw)
##no dispersal
wsimp1nd<-glmmTMB(simp.norm~day+temperature+infusion+(1|mesocosm), data=df.rare.div.mw)
##minus infusion
wsimp1ni<-glmmTMB(simp.norm~day+temperature+dispersal+(1|mesocosm), data=df.rare.div.mw)
##minus temperature
wsimp1nt<-glmmTMB(simp.norm~day+infusion+dispersal+(1|mesocosm), data=df.rare.div.mw)
##minus time
wsimp1nnd<-glmmTMB(simp.norm~temperature+infusion+dispersal+(1|mesocosm), data=df.rare.div.mw)
anova(wsimp1,wsimp1nd) #0.87## Data: df.rare.div.mw
## Models:
## wsimp1nd: simp.norm ~ day + temperature + infusion + (1 | mesocosm), zi=~0, disp=~1
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wsimp1nd 8 385.60 412.41 -184.80 369.60
## wsimp1 9 387.57 417.74 -184.79 369.57 0.0268 1 0.8699## Data: df.rare.div.mw
## Models:
## wsimp1ni: simp.norm ~ day + temperature + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wsimp1ni 7 487.79 511.25 -236.89 473.79
## wsimp1 9 387.57 417.74 -184.79 369.57 104.21 2 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: df.rare.div.mw
## Models:
## wsimp1nt: simp.norm ~ day + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wsimp1nt 8 418.72 445.53 -201.36 402.72
## wsimp1 9 387.57 417.74 -184.79 369.57 33.144 1 8.56e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: df.rare.div.mw
## Models:
## wsimp1nnd: simp.norm ~ temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wsimp1nnd 7 427.04 450.50 -206.52 413.04
## wsimp1 9 387.57 417.74 -184.79 369.57 43.464 2 3.646e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: simp.norm
## Chisq Df Pr(>Chisq)
## day 50.7854 2 9.377e-12 ***
## temperature 42.1861 1 8.299e-11 ***
## infusion 231.0011 2 < 2.2e-16 ***
## dispersal 0.0268 1 0.87
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Analysis of Deviance Table (Type II Wald chisquare tests)
#
# Response: simp.norm
# Chisq Df Pr(>Chisq)
# day 50.7856 2 9.377e-12 ***
# temperature 42.1866 1 8.297e-11 ***
# infusion 231.0025 2 < 2.2e-16 ***
# dispersal 0.0268 1 0.87
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##interactions?
wsimp1.int <- glmmTMB(simp.norm~day*infusion+temperature+dispersal+(1|mesocosm),data=df.rare.div.mw)
wsimp1.int1 <- glmmTMB(simp.norm~day+infusion*temperature+dispersal+(1|mesocosm),data=df.rare.div.mw)
wsimp1.int2 <- glmmTMB(simp.norm~day+infusion+temperature*dispersal+(1|mesocosm),data=df.rare.div.mw)
wsimp1.int3 <- glmmTMB(simp.norm~day*temperature+infusion+dispersal+(1|mesocosm),data=df.rare.div.mw)
wsimp1.int4 <- glmmTMB(simp.norm~day*dispersal+temperature+infusion+(1|mesocosm),data=df.rare.div.mw)
wsimp1.int5 <- glmmTMB(simp.norm~day+dispersal*infusion+temperature+(1|mesocosm),data=df.rare.div.mw)
anova(wsimp1,wsimp1.int) #1.148e-09***## Data: df.rare.div.mw
## Models:
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1.int: simp.norm ~ day * infusion + temperature + dispersal + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wsimp1 9 387.57 417.74 -184.79 369.57
## wsimp1.int 13 347.98 391.56 -160.99 321.98 47.592 4 1.148e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: df.rare.div.mw
## Models:
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1.int1: simp.norm ~ day + infusion * temperature + dispersal + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wsimp1 9 387.57 417.74 -184.79 369.57
## wsimp1.int1 11 387.64 424.51 -182.82 365.64 3.9376 2 0.1396## Data: df.rare.div.mw
## Models:
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1.int2: simp.norm ~ day + infusion + temperature * dispersal + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wsimp1 9 387.57 417.74 -184.79 369.57
## wsimp1.int2 10 388.86 422.38 -184.43 368.86 0.7141 1 0.3981## Data: df.rare.div.mw
## Models:
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1.int3: simp.norm ~ day * temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wsimp1 9 387.57 417.74 -184.79 369.57
## wsimp1.int3 11 391.02 427.89 -184.51 369.02 0.5552 2 0.7576## Data: df.rare.div.mw
## Models:
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1.int4: simp.norm ~ day * dispersal + temperature + infusion + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wsimp1 9 387.57 417.74 -184.79 369.57
## wsimp1.int4 11 389.79 426.66 -183.90 367.79 1.7812 2 0.4104## Data: df.rare.div.mw
## Models:
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1.int5: simp.norm ~ day + dispersal * infusion + temperature + (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wsimp1 9 387.57 417.74 -184.79 369.57
## wsimp1.int5 11 391.00 427.87 -184.50 369.00 0.5773 2 0.7493## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: simp.norm
## Chisq Df Pr(>Chisq)
## day 71.6034 2 2.828e-16 ***
## infusion 229.6298 2 < 2.2e-16 ***
## temperature 40.5261 1 1.940e-10 ***
## dispersal 0.0345 1 0.8526
## day:infusion 57.0587 4 1.203e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Analysis of Deviance Table (Type II Wald chisquare tests)
#
# Response: simp.norm
# Chisq Df Pr(>Chisq)
# day 71.6032 2 2.829e-16 ***
# infusion 229.6270 2 < 2.2e-16 ***
# temperature 40.5264 1 1.940e-10 ***
# dispersal 0.0345 1 0.8526
# day:infusion 57.0591 4 1.202e-11 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#tapply(df.rare.div.mw$log(simp.norm), df.rare.div.mw$infusion, mean)
#tapply(df.rare.div.mw$log(simp.norm), df.rare.div.mw$day, mean)
##model checking
shapiro.test(residuals(wsimp1.int)) #awesome##
## Shapiro-Wilk normality test
##
## data: residuals(wsimp1.int)
## W = 0.93335, p-value = 3.21e-08
##posthoc things
##just without dispersal because not plotting it
##also not plotting day but it should be in the model somewhere already
wsimp1.int.em <- emmeans(wsimp1.int,~infusion+temperature)
multcomp::cld(wsimp1.int.em)## infusion temperature emmean SE df lower.CL upper.CL .group
## SG C -1.2646 0.0937 198 -1.449 -1.0798 1
## SG H -0.6694 0.0941 198 -0.855 -0.4839 2
## OL C -0.0993 0.0931 198 -0.283 0.0843 3
## SW C 0.4376 0.0931 198 0.254 0.6213 4
## OL H 0.4959 0.0935 198 0.312 0.6803 4
## SW H 1.0328 0.0935 198 0.848 1.2172 5
##
## Results are averaged over the levels of: day, dispersal
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 6 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.3.3.2.3 Mosquito Simpson’s
From manuscript pre-revisions: “In contrast, the Simpson’s index, which integrates evenness into a measure of alpha diversity, of the adult mosquito microbiome was lower in male mosquitoes relative to female mosquitoes (p<0.0001). This trend was associated with the increased dominance of wAlbB in males. The Simpson’s index of the adult mosquito microbiome did not vary significantly with the dispersal, aquatic chemistry, temperature, or time of emergence.”
##full model
msimp.all <- glmmTMB(simp.norm~newday+temperature+infusion+dispersal+sex+(1|mesocosm), data=df.div.mq)
##without sex
msimp.nosex <- glmmTMB(simp.norm~newday+temperature+infusion+dispersal+(1|mesocosm), data=df.div.mq)
##without dispersal
msimp.nodis <- glmmTMB(simp.norm~newday+temperature+infusion+sex+(1|mesocosm), data=df.div.mq)
##without infusion
msimp.noinf <- glmmTMB(simp.norm~newday+temperature+dispersal+sex+(1|mesocosm), data=df.div.mq)
##without temperature
msimp.notem <- glmmTMB(simp.norm~newday+infusion+dispersal+sex+(1|mesocosm), data=df.div.mq)
##without time
msimp.notim <- glmmTMB(simp.norm~temperature+infusion+dispersal+sex+(1|mesocosm), data=df.div.mq)
anova(msimp.all,msimp.nodis) #0.10## Data: df.div.mq
## Models:
## msimp.nodis: simp.norm ~ newday + temperature + infusion + sex + (1 | mesocosm), zi=~0, disp=~1
## msimp.all: simp.norm ~ newday + temperature + infusion + dispersal + sex + , zi=~0, disp=~1
## msimp.all: (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## msimp.nodis 9 507.26 536.72 -244.63 489.26
## msimp.all 10 506.57 539.30 -243.28 486.57 2.6991 1 0.1004## Data: df.div.mq
## Models:
## msimp.noinf: simp.norm ~ newday + temperature + dispersal + sex + (1 | mesocosm), zi=~0, disp=~1
## msimp.all: simp.norm ~ newday + temperature + infusion + dispersal + sex + , zi=~0, disp=~1
## msimp.all: (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## msimp.noinf 8 505.29 531.47 -244.65 489.29
## msimp.all 10 506.57 539.30 -243.28 486.57 2.7247 2 0.2561## Data: df.div.mq
## Models:
## msimp.nosex: simp.norm ~ newday + temperature + infusion + dispersal + (1 | , zi=~0, disp=~1
## msimp.nosex: mesocosm), zi=~0, disp=~1
## msimp.all: simp.norm ~ newday + temperature + infusion + dispersal + sex + , zi=~0, disp=~1
## msimp.all: (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## msimp.nosex 9 562.49 591.95 -272.25 544.49
## msimp.all 10 506.57 539.30 -243.28 486.57 57.925 1 2.723e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: df.div.mq
## Models:
## msimp.notem: simp.norm ~ newday + infusion + dispersal + sex + (1 | mesocosm), zi=~0, disp=~1
## msimp.all: simp.norm ~ newday + temperature + infusion + dispersal + sex + , zi=~0, disp=~1
## msimp.all: (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## msimp.notem 9 504.87 534.33 -243.44 486.87
## msimp.all 10 506.57 539.30 -243.28 486.57 0.306 1 0.5801## Data: df.div.mq
## Models:
## msimp.notim: simp.norm ~ temperature + infusion + dispersal + sex + (1 | mesocosm), zi=~0, disp=~1
## msimp.all: simp.norm ~ newday + temperature + infusion + dispersal + sex + , zi=~0, disp=~1
## msimp.all: (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## msimp.notim 8 504.93 531.12 -244.47 488.93
## msimp.all 10 506.57 539.30 -243.28 486.57 2.3679 2 0.3061## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: simp.norm
## Chisq Df Pr(>Chisq)
## newday 2.3823 2 0.30387
## temperature 0.3062 1 0.58000
## infusion 2.7438 2 0.25363
## dispersal 2.8697 1 0.09026 .
## sex 67.4478 1 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Response: simp.norm
# Chisq Df Pr(>Chisq)
# newday 2.3823 2 0.30387
# temperature 0.3062 1 0.58000
# infusion 2.7438 2 0.25363
# dispersal 2.8697 1 0.09026 .
# sex 67.4477 1 < 2e-16 ***
#AICtab(msimp.all,msimp.nosex,msimp.nodis,msimp.noinf,msimp.notem,msimp.notim)
#no time & no temp tied for best, then no infusion
##all of the below was bad before normalizing
shapiro.test(residuals(msimp.all))##
## Shapiro-Wilk normality test
##
## data: residuals(msimp.all)
## W = 0.98791, p-value = 0.0961
##posthoc things
msimp.all.em <- emmeans(msimp.all,~infusion*temperature)
multcomp::cld(msimp.all.em)## infusion temperature emmean SE df lower.CL upper.CL .group
## SG H -0.2061 0.152 185 -0.506 0.0939 1
## SG C -0.1150 0.167 185 -0.445 0.2150 1
## OL H -0.0782 0.143 185 -0.361 0.2045 1
## OL C 0.0129 0.115 185 -0.215 0.2406 1
## SW H 0.1484 0.162 185 -0.172 0.4682 1
## SW C 0.2394 0.193 185 -0.141 0.6197 1
##
## Results are averaged over the levels of: newday, dispersal, sex
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 6 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.3.3.2.4 Mosquito Simpson’s rarefied
df.rare.div.mq$simp.norm <- bestNormalize(df.rare.div.mq$InvSimpson)$x.t #orderNorm
hist(df.rare.div.mq$simp.norm)
##full model
msimp.all <- glmmTMB(simp.norm~newday+temperature+infusion+dispersal+sex+(1|mesocosm), data=df.rare.div.mq)
##without sex
msimp.nosex <- glmmTMB(simp.norm~newday+temperature+infusion+dispersal+(1|mesocosm), data=df.rare.div.mq)
##without dispersal
msimp.nodis <- glmmTMB(simp.norm~newday+temperature+infusion+sex+(1|mesocosm), data=df.rare.div.mq)
##without infusion
msimp.noinf <- glmmTMB(simp.norm~newday+temperature+dispersal+sex+(1|mesocosm), data=df.rare.div.mq)
##without temperature
msimp.notem <- glmmTMB(simp.norm~newday+infusion+dispersal+sex+(1|mesocosm), data=df.rare.div.mq)
##without time
msimp.notim <- glmmTMB(simp.norm~temperature+infusion+dispersal+sex+(1|mesocosm), data=df.rare.div.mq)
anova(msimp.all,msimp.nodis) #0.095 .## Data: df.rare.div.mq
## Models:
## msimp.nodis: simp.norm ~ newday + temperature + infusion + sex + (1 | mesocosm), zi=~0, disp=~1
## msimp.all: simp.norm ~ newday + temperature + infusion + dispersal + sex + , zi=~0, disp=~1
## msimp.all: (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## msimp.nodis 9 506.90 536.36 -244.45 488.90
## msimp.all 10 506.12 538.85 -243.06 486.12 2.7861 1 0.09508 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: df.rare.div.mq
## Models:
## msimp.noinf: simp.norm ~ newday + temperature + dispersal + sex + (1 | mesocosm), zi=~0, disp=~1
## msimp.all: simp.norm ~ newday + temperature + infusion + dispersal + sex + , zi=~0, disp=~1
## msimp.all: (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## msimp.noinf 8 504.89 531.07 -244.44 488.89
## msimp.all 10 506.12 538.85 -243.06 486.12 2.7719 2 0.2501## Data: df.rare.div.mq
## Models:
## msimp.nosex: simp.norm ~ newday + temperature + infusion + dispersal + (1 | , zi=~0, disp=~1
## msimp.nosex: mesocosm), zi=~0, disp=~1
## msimp.all: simp.norm ~ newday + temperature + infusion + dispersal + sex + , zi=~0, disp=~1
## msimp.all: (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## msimp.nosex 9 562.45 591.91 -272.23 544.45
## msimp.all 10 506.12 538.85 -243.06 486.12 58.334 1 2.212e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: df.rare.div.mq
## Models:
## msimp.notem: simp.norm ~ newday + infusion + dispersal + sex + (1 | mesocosm), zi=~0, disp=~1
## msimp.all: simp.norm ~ newday + temperature + infusion + dispersal + sex + , zi=~0, disp=~1
## msimp.all: (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## msimp.notem 9 504.45 533.91 -243.22 486.45
## msimp.all 10 506.12 538.85 -243.06 486.12 0.333 1 0.5639## Data: df.rare.div.mq
## Models:
## msimp.notim: simp.norm ~ temperature + infusion + dispersal + sex + (1 | mesocosm), zi=~0, disp=~1
## msimp.all: simp.norm ~ newday + temperature + infusion + dispersal + sex + , zi=~0, disp=~1
## msimp.all: (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## msimp.notim 8 504.45 530.64 -244.23 488.45
## msimp.all 10 506.12 538.85 -243.06 486.12 2.3387 2 0.3106## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: simp.norm
## Chisq Df Pr(>Chisq)
## newday 2.3527 2 0.30840
## temperature 0.3332 1 0.56376
## infusion 2.7917 2 0.24762
## dispersal 2.9784 1 0.08438 .
## sex 67.9982 1 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Chisq Df Pr(>Chisq)
# newday 2.3527 2 0.30840
# temperature 0.3332 1 0.56376
# infusion 2.7917 2 0.24762
# dispersal 2.9784 1 0.08438 .
# sex 67.9982 1 < 2e-16 ***
#AICtab(msimp.all,msimp.nosex,msimp.nodis,msimp.noinf,msimp.notem,msimp.notim)
#no time & no temp tied for best, then no infusion
##all of the below was bad before normalizing
shapiro.test(residuals(msimp.all))##
## Shapiro-Wilk normality test
##
## data: residuals(msimp.all)
## W = 0.98823, p-value = 0.107
##posthoc things
msimp.all.em <- emmeans(msimp.all,~infusion*temperature)
multcomp::cld(msimp.all.em)## infusion temperature emmean SE df lower.CL upper.CL .group
## SG H -0.2117 0.152 185 -0.511 0.0879 1
## SG C -0.1168 0.167 185 -0.446 0.2128 1
## OL H -0.0779 0.143 185 -0.360 0.2045 1
## OL C 0.0170 0.115 185 -0.210 0.2445 1
## SW H 0.1446 0.162 185 -0.175 0.4641 1
## SW C 0.2395 0.193 185 -0.140 0.6193 1
##
## Results are averaged over the levels of: newday, dispersal, sex
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 6 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.3.4 Infusion water things
df.div.inf <- subset(df.div,type=="Infusion water")
df.div.inf.se <- summarySE(df.div.inf,measurevar="Observed",groupvars=c("infusion"))
df.div.inf.se## infusion N Observed sd se ci
## 1 OL 3 45.666667 3.214550 1.8559215 7.985386
## 2 SG 3 23.666667 6.350853 3.6666667 15.776393
## 3 SW 3 4.333333 1.154701 0.6666667 2.868435# gg.iw.rich <- ggplot(df.div.inf,aes(x=infusion,y=Observed,fill=infusion))+
# geom_jitter(alpha=0.5,position=position_jitterdodge(jitter.width=0.2),color="gray")+
# geom_errorbar(data=df.div.inf.se,aes(ymax=Observed+se,ymin=Observed-se),width=0.2,color="black",position=position_dodge(width=0.6))+
# geom_point(data=df.div.inf.se,size=2.5,position=position_dodge(width=0.6))+
# #geom_boxplot(outlier.shape=NA,alpha=0.5)+
# theme_cowplot()+
# ggtitle("Infusion water")
# gg.iw.rich4 Diversity correlations
4.1 Diversity x day
Doesn’t look interesting
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
4.2 Diversity - mesocosm vs. mosquito
##
## Attaching package: 'reshape'## The following object is masked from 'package:cowplot':
##
## stamp## The following objects are masked from 'package:plyr':
##
## rename, round_anymeso.corr <- merge(df.div.mq,df.div.mw,by="mesocosm")
ggplot(meso.corr,aes(x=InvSimpson.x,y=InvSimpson.y))+
geom_point()+
#facet_wrap(~day.y)+
geom_smooth()## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
5 Session info
## R version 4.4.1 (2024-06-14)
## Platform: aarch64-apple-darwin20
## Running under: macOS Sonoma 14.2.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: Pacific/Honolulu
## tzcode source: internal
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] reshape_0.8.9 bestNormalize_1.9.1 DHARMa_0.4.6
## [4] bbmle_1.0.25.1 emmeans_1.10.3 glmmTMB_1.1.9
## [7] phyloseq_1.48.0 cowplot_1.1.3 ggpubr_0.6.0
## [10] ggplot2_3.5.1 Rmisc_1.5.1 plyr_1.8.9
## [13] lattice_0.22-6 car_3.1-2 carData_3.0-5
##
## loaded via a namespace (and not attached):
## [1] jsonlite_1.8.8 magrittr_2.0.3 TH.data_1.1-2
## [4] estimability_1.5.1 farver_2.1.2 nloptr_2.1.1
## [7] rmarkdown_2.27 zlibbioc_1.50.0 vctrs_0.6.5
## [10] multtest_2.60.0 minqa_1.2.7 rstatix_0.7.2
## [13] butcher_0.3.4 htmltools_0.5.8.1 broom_1.0.6
## [16] Rhdf5lib_1.26.0 rhdf5_2.48.0 sass_0.4.9
## [19] parallelly_1.38.0 bslib_0.8.0 sandwich_3.1-0
## [22] lubridate_1.9.3 zoo_1.8-12 cachem_1.1.0
## [25] TMB_1.9.14 igraph_2.0.3 mime_0.12
## [28] lifecycle_1.0.4 iterators_1.0.14 pkgconfig_2.0.3
## [31] gap_1.5-3 Matrix_1.7-0 R6_2.5.1
## [34] fastmap_1.2.0 rbibutils_2.2.16 shiny_1.9.1
## [37] GenomeInfoDbData_1.2.12 future_1.34.0 digest_0.6.36
## [40] numDeriv_2016.8-1.1 colorspace_2.1-1 S4Vectors_0.42.1
## [43] qgam_1.3.4 vegan_2.6-6.1 labeling_0.4.3
## [46] timechange_0.3.0 fansi_1.0.6 httr_1.4.7
## [49] abind_1.4-5 mgcv_1.9-1 compiler_4.4.1
## [52] rngtools_1.5.2 withr_3.0.1 doParallel_1.0.17
## [55] backports_1.5.0 highr_0.11 ggsignif_0.6.4
## [58] MASS_7.3-61 lava_1.8.0 biomformat_1.32.0
## [61] permute_0.9-7 tools_4.4.1 ape_5.8
## [64] httpuv_1.6.15 future.apply_1.11.2 nnet_7.3-19
## [67] glue_1.7.0 promises_1.3.0 nlme_3.1-165
## [70] rhdf5filters_1.16.0 grid_4.4.1 cluster_2.1.6
## [73] reshape2_1.4.4 ade4_1.7-22 generics_0.1.3
## [76] recipes_1.1.0 gtable_0.3.5 nortest_1.0-4
## [79] class_7.3-22 tidyr_1.3.1 data.table_1.15.4
## [82] utf8_1.2.4 XVector_0.44.0 BiocGenerics_0.50.0
## [85] foreach_1.5.2 pillar_1.9.0 stringr_1.5.1
## [88] later_1.3.2 splines_4.4.1 dplyr_1.1.4
## [91] survival_3.7-0 tidyselect_1.2.1 Biostrings_2.72.1
## [94] knitr_1.48 gridExtra_2.3 IRanges_2.38.1
## [97] xfun_0.46 Biobase_2.64.0 hardhat_1.4.0
## [100] timeDate_4032.109 stringi_1.8.4 UCSC.utils_1.0.0
## [103] yaml_2.3.10 boot_1.3-30 evaluate_0.24.0
## [106] codetools_0.2-20 tibble_3.2.1 multcompView_0.1-10
## [109] cli_3.6.3 rpart_4.1.23 Rdpack_2.6
## [112] xtable_1.8-4 munsell_0.5.1 jquerylib_0.1.4
## [115] Rcpp_1.0.13 GenomeInfoDb_1.40.1 globals_0.16.3
## [118] bdsmatrix_1.3-7 parallel_4.4.1 gower_1.0.1
## [121] gap.datasets_0.0.6 doRNG_1.8.6 lme4_1.1-35.5
## [124] listenv_0.9.1 mvtnorm_1.2-5 ipred_0.9-15
## [127] scales_1.3.0 prodlim_2024.06.25 purrr_1.0.2
## [130] crayon_1.5.3 rlang_1.1.4 multcomp_1.4-26